Matter- wave solitons in the counterflow of two immiscible superfluids 
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We study formation of solitons induced by counterflows of immiscible superfluids. Our setting is based on a 
quasi-one-dimensional binary Bose-Einstein condensate (BEC), composed of two immiscible components with 
large and small numbers of atoms in them. Assuming that the "small" component moves with constant velocity, 
either by itself, or being dragged by a moving trap, and intrudes into the "large" counterpart, the following 
results are obtained. Depending on the velocity, and on whether the small component moves in the absence or in 
the presence of the trap, two-component dark-bright solitons, scalar dark solitons, or multiple dark solitons may 
emerge, the latter outcome taking place due to breakdown of the superfluidity. We present two sets of analytical 
results to describe this phenomenology. In an intermediate velocity regime, where dark-bright solitons form, 
a reduction of the two-component Gross-Pitaevskii system to an integrable Mel'nikov system is developed, 
demonstrating that solitary waves of the former are very accurately described by analytically available solitons 
of the latter. In the high- velocity regime, where the breakdown of the superfluidity induces the formation of 
dark solitons and multi-soliton trains, an effective single-component description, in which a strongly localized 
wave packet of the "small" component acts as an effective potential for the "large" one, allows us to estimate 
the critical velocity beyond which the coherent structures emerge in good agreement with the numerical results. 
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I. INTRODUCTION 



Matter-wave solitons is a theme that has attracted much 
attention in studies of atomic Bose-Einstein condensates 
(BECs) d H and, more generally, in studies of superflu- 
ids 0] - see, e.g., recent reviews In atomic BECs, 
matter- wave solitons have been realized by means of various 
methods, including phase-imprinting [7], density engineering 
[80, quantum- state engineering 09D, transport in optical lattices 
I Toll , transport past single defects [Hill or through disordered 
potentials II 1 211 . and nonlinear interference between individ- 
ual condensates 11311 . Furthermore, matter- wave solitons have 
also been observed in multi-component BECs consisting of 
two different spin states of the same atom species. In partic- 
ular, atomic dark-bright solitons have been realized in binary 
rubidium condensates by means of the quantum- state engi- 
neering technique |9D, or the counterflow between two misci- 
ble components [lMLL5l]; employing the latter technique, gen- 
eration of dark-dark solitons was reported too Ql^l Ol • Im- 
portantly, solitons have been argued to emerge spontaneously 
upon crossing the BEC phase transition, via the Kibble-Zurek 
mechanism 1181 . This can be thought of as a quasi-one- 
dimensional analog of the pioneering experiments on the same 
mechanism in higher-dimensional BECs performed in [fl9ll . 

As mentioned above, the counterflow of two miscible, 
quasi-one-dimensional (ID) BEC superfluids, is one of the 
techniques that were recently used for the generation of 
matter- wave solitons. More generally, the counterflow dy- 
namics of superfluids is an interesting subject due to its links 
with such fundamental phenomena as turbulence, pattern for- 
mation, and generation of topological defects, which occur 
in a vast variety of physical contexts ranging from superfluid 
helium to cosmological strings and inflation [20] . In the con- 
text of binary atomic BECs, a number of theoretical works 



have been dealing with the counterflow of either miscible 
IH [23] or immiscible HQ binary superfluids. In the for- 
mer case, quantum turbulence and vortex generation in higher- 
dimensional settings i2lll and in spinor BECs with F = 1 i22ll 
were studied; in the latter case, effects of quantum swapping, 
Rayleigh-Taylor instability l23ll . and parametric resonance of 
capillary waves were predicted [24]. Furthermore, more re- 
cently, the counterflow superfluid of polaron pairs in Bose- 
Fermi mixtures in optical lattices was analyzed [|25ll . 

In the present work, we propose and analyze a physically 
relevant setting based on the counterflow of two immiscible 
superfluids. This setting can be used for the generation of 
scalar or vector solitons (of the dark-bright type) in binary 
BECs, and also for studies of their superfluid properties. We 
consider a quasi- ID binary BEC mixture composed of two 
immiscible components, "small" and "large" ones, with ratios 
of axial sizes and number of atoms in them ~ 1 : 100. We 
then assume that the "small" component moves with constant 
velocity v through the "large" one, by itself or being dragged 
by a moving trapping potential. 

Systematic simulations of this setting produce the follow- 
ing phenomenology. For sufficiently small dragging veloci- 
ties, v < 0.61c s , with c s = V2c s being the speed of sound of 
the large component, and if the trap of the small component 
is switched off, the small component cannot remain localized, 
hence no well-defined patterns are formed. For intermediate 
velocities, 0.61c s < v < 0.95c s , a robust two-component 
dark-bright soliton is formed. It is composed of a bright (dark) 
soliton in the small (large) component and can be very well 
approximated analytically. In particular, we use a multiscale 
asymptotic expansion method to reduce the underlying system 
of the Gross-Pitaevskii equations (GPEs) to a completely inte- 
grable Mel'nikov system [12611. thus very accurately mapping 
the dark-bright solitary wave, which spontaneously emerges 
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in the GPE system, into an exact soliton of the Mel'nikov 
system. On the other hand, for 0.58c s < v < 0.87c s but 
in the presence of its trap dragging the small component, the 
"Mel'nikov's" dark-bright soliton still emerges but is unstable 
and eventually decays. For larger trap velocities, v > 0.95c s , 
and if the trap is absent, a scalar dark soliton is created in the 
large component. Furthermore, for v > 0.87c s and in the 
presence of the trap, the small component is deformed into 
a sharply localized object, with a width on the order of the 
healing length, which propagates with a velocity larger than 
the critical velocity for the breakdown of the superfluidity in 
the large component. The breakdown leads to the emergence 
of multiple dark solitons, similarly to the results reported in 
previous theoretical J27|,|28|] and experimental [11] works (cf. 
also the quasi-two-dimensional analog of such an experiment 
producing vortex dipoles fl29ID . 

The rest of the paper is structured as follows. The model 
is formulated in Section II. Then, in Section III, we system- 
atically analyze, by means of analytical approximations and 
direct simulations, two different scenarios of the evolution of 
the small component, in the absence or in the presence of its 
trap. Finally, in Section IV we summarize our findings and 
discuss possibilities for the extension of the analysis. 

II. THE MODEL AND SETUP 

We consider the binary BEC composed of two hyperfine 
states of 87 Rb, that are described by macroscopic wave func- 
tions iffj (r, t) = 1, 2). The dynamics can then be described 
by the following system of coupled GPEs 01 0] : 

ihdt^j = (^-^V 2 + Vj(r) - fjLj + X>i*|tf*| 2 ^ *i> 

(1) 

where m is the atomic mass, fij are chemical potentials, Vj (r) 
are trapping potentials, and gjk = 47rft 2 a 7 / c /m are coupling 
constants defined by the s-wave scattering lengths, a^. No- 
tice that, for states |1, — 1) and |2, 1), which were used in the 
experiment of Ref. fl30ll . or states |1, — 1) and |2, — 2) used in 
Refs. lfl~4L ITsh , all the scattering lengths take approximately 
equal (positive) values. Below we assume an = a 22 7^ 
ai2, hence the sign of parameter A = (an 022 — a i2)/ a ii' 
which determines whether the two components are miscible 
(A > 0) or immiscible (A < 0) [31], depends only on ratio 
f3 = a\ 2 ja\\\ if f3 > 1 (/? < 1), then the two superfluids are 
immiscible (miscible). 

We consider the case when both components are confined 
in strongly anisotropic (quasi- ID) traps, which have the form 
of rectangular boxes of lengths L Xj ^> L Vj = L Zj = L±, 
with the transverse-confinement lengths L± j being on the or- 
der of the healing lengths £j . For such a highly anisotropic 
trap, it is relevant to factorize the wave functions into longi- 
tudinal and transverse components, namely {u(x, i),v(x, t)} 
and $j(y, z) (see, e.g., Refs. 0-11]): 

*i(r,t) = $i(2/,2?Xa;,t)exp(-zJ5;it/fi), (2) 
*2(r,t) = <S> 2 (y,z)v(x : t)exp(-iE 2 t/h), (3) 
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FIG. 1. (Color online) The density profiles (solid lines) of the two 
immiscible superfluids, confined by the box-shaped traps V± (x) and 
Vi{x, t) (dashed lines) at t = 0. The red (green) line corresponds to 
the u- (v-) component. The arrow depicts the direction of motion of 
the v-component at t > 0. 

where the transverse modal functions <$>j(y,z) and ener- 
gies Ej are determined by the auxiliary 2D problem for 
the quantum oscillator in the transverse box. Since the 
considered system is effectively one-dimensional, it is nat- 
ural to assume that &j remains in the ground state, i.e., 
®j(y,z) = (V2/L±)sm(7ry/L_ L )sm(7rz/L±) 9 forj = 1,2. 
Using these expressions, we substitute expressions ©-(13]) into 
Eqs. (Q]) and perform averaging over the transverse coordi- 
nates, to derive a system of scaled equations for the longi- 
tudinal wave functions u(x, t) and v(x, t) (see also Ref. 13211 ): 

iu t = -u xx + (M 2 + /3|v| 2 - /ii> + Vi(x)u, (4) 
ivt = -v xx + WH 2 + |v| 2 - /i 2 )v + V 2 (x)v, (5) 

where subscripts stand for partial derivatives. In these equa- 
tions, longitudinal coordinate x, time t, densities \u\ 2 , |v| 2 , 
and energy are measured, respectively, in units of the heal- 
ing length £1 = h/ y/2mn\gii, characteristic time V^^i/cs 
(where c s = \Jg\\n\lm is the sound speed of the first com- 
ponent), density n\, and characteristic energy g\\ri\ of the 
first component. Finally, g^ = (9/41^ )g^ are the effective 
ID interaction strengths. 

It is now useful to adopt experimentally relevant values of 
the parameters. First, we fix the scattering length and the den- 
sity of first component to an = 5.77 nm and m = 10 8 m _1 , 
which results in the healing length, £1 = 0.25 /mi; next, we 
fix f3 = 1.1, which corresponds to the immiscibility (we have 
checked that other values of f3 > 1 lead to qualitatively similar 
results). We will assume that the first component (u) is "much 
larger" than the other one (v), in the sense that the respective 
ratio of the numbers of atoms is Ni/N 2 ~ 100; accordingly, 
the chemical potentials are fixed to [i\ = 1 and [i 2 = 0.06. 
As concerns the size of the traps, we assume that, in most 
cases, L Xl = 57.6 //m, L X2 = 5 //m, and L± = 0.8 /im, 
which correspond, respectively, to L Xl = 244, L X2 =21 and 
L± = 3 in the dimensionless units; an exception concerns 
Fig. [5] (see below), where we use L Xl = 123 /im (or 490 in 
the dimensionless units) for the length of the li-component. 
We have checked that, as long as the ratios of the atom 
numbers, chemical potentials and trap lengths are such that 
N 1 /N 2 - /ii//i 2 - 100 and L Xl : L X2 : L± - 100 : 10 : 1, 
the results are qualitatively similar. 
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FIG. 2. (Color online) Contour plots showing the evolution of 
the largest density dip in the ^-component (the top panel) and the 
peak density of the v-component (the bottom panel), with respect 
to the velocity of the motion of the trap. We identify three dif- 
ferent regimes, separated by the vertical (green) dashed lines: a) 
< v < 0.61, where no localized excitations are formed, b) 
0.61 < v < 0.95, where traveling dark-bright solitons emerge, 
and c) v > 0.95, where deep dark solitons are formed in the u- 
component. In the latter case, the v-component is quickly dispersed. 



In our simulations, we will approximate the box-like trap- 
ping potentials Vj by super-Gaussians: 



Vj = V 



1 — exp 



X + Xj 



Wj 



12\ 



(6) 



where Vo is the (common) trap amplitude for both compo- 
nents, while Xj and Wj denote the positions and widths of the 
traps, respectively. We fix these parameters as Vq = 10, x\ = 
—20, w\ = 160, and W2 = 20; as concerns the position of the 
trap acting on the v-component, it is assumed to be initially (at 
t = 0) placed at X2 (0) = 135 and move at a constant velocity: 
X2 — 135 — vt. Thus, the "small" v-component starts moving 
to the right, penetrating the "large" ^-component; notice that 
v is measured in units of c s — c s /V2 (recall that c s is the 
speed of sound of the ^-component). As we show below, such 
an effective counterflow gives rise to the formation of nonlin- 
ear structures and solitons that we will investigate in detail in 
the following section. Furthermore, we note in passing that 
in the immiscible case considered here, the counterflow with 
speed v is modulationally stable, in accordance with the cri- 
terion of EH EH]- This can also be intuitively understood on 
the basis of immiscibility of two independently modulation- 
ally stable components. 

The initial configuration, obtained by means of an 
imaginary-time integration, is depicted in Fig. [T] where the 



two components are initially separated. Then, at t > 0, the 
v-component is set into motion and moves through the u- 
component. The evolution of the system is monitored by nu- 
merically integrating Eqs. © and © in real time by means of 
the split-step Fourier method. 



III. ANALYTICAL AND NUMERICAL RESULTS 

We will consider two different cases: (A) when trap V2, act- 
ing on the small component, is switched off after it has pene- 
trated the ^-component (for the above-mentioned parameters, 
this occurs at x = —40 and t ~ 100); (B) the trap V2 is not 
switched off, i.e., the v-component is permanently dragged by 
its trap. 



A. Evolution of the small component in the absence of the trap. 

In this case, our systematic simulations have revealed the 
existence of three different velocity regimes, characterized by 
the emergence (or the failure to emerge) of various localized 
nonlinear structures. These regimes were identified by calcu- 
lating - at each time step - the maximum density dip in the 
^-component, as well as the peak density of the v-component, 
as shown in the top and bottom panel of Fig. O respectively. 

In particular, for velocities < v < 0.61 [correspond- 
ing to the region on the left of the first dashed (green) line 
in Fig. O, no clear solitary- wave structure is formed in the 
dynamics. There exist, however, numerous small- amplitude 
dips that are formed in the ^-component, which propagate 
with speeds approximately equal to the speed of sound in that 
component. Generally, for such relatively small velocities, as 
long as the v-component penetrates the ^-component, strong 
emission of radiation is observed. A typical example is shown 
in Fig. [3j for v = 0.3. In addition to the small- amplitude 
sound wavepackets propagating in the u component, there ex- 
ists a slow structure (with a speed near close to that of the 
initially moving trap), which emerges with a larger dip in the 
^-component and a bright "bump" in the v component. Yet, 
this structure does not preserve its shape, but rather slowly dis- 
perses, contrary to what would be expected of a robust solitary 
wave. 

On the other hand, in the intermediate velocity regime, 
0.61 < v < 0.95 [see the region between dashed (green) 
lines in Fig. [2), the v-component forms a relatively large peak, 
which is practically robust and does not decay in time. Addi- 
tionally, a dip of a constant density is formed in the large (u) 
component, which propagates with the same velocity as the 
peak in the the v-component, i.e., the two components build 
a traveling coherent structure. An example, shown in the top 
panel of Fig. [4] for v — 0.7, reveals this in more detail: after 
a relatively short transient period (t ~ 100) needed for the 
v-component to reach the flat segment of the ^-component, a 
small-amplitude dark-bright soliton is formed. Notice that the 
strong emission of sound waves -which is observed chiefly in 
the ^-component- is solely due to the penetration process (as 
in the case of Fig. 0: the sound waves get detached from the 



FIG. 3. (Color online) Contour plots showing the evolution of the 
density of the ^-component (the inset shows the evolution of the v- 
component) for trap's velocity v = 0.3. It is observed, aside from 
the emission of sound waves in the ^-component, that a structure 
resembling a dark-bright solitary wave is formed. However, it slowly 
disperses in the course of the evolution. 



FIG. 5. (Color online) The same as in the middle and bottom panels 
of Fig. HI but for a larger size of trap Vi, namely L X1 = 123 pm 
(or L Xl — 490 in dimensionless units). In this case, as shown in the 
snapshots at t = 200 (the top panel) and t = 376 (the bottom panel), 
the evolution time of the Mel'nikov soliton is significantly longer. 
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FIG. 4. (Color online) The top panel shows the evolution of the den- 
sity the ^-component (inset shows the evolution of the v-component) 
for trap's velocity v = 0.7. It is observed that a Mel'nikov's soliton, 
moving with a velocity ^ so i = 0.96, is formed. The middle and bot- 
tom panels show snapshots of the densities of the dark soliton in the 
^-component, and the bright soliton in the v-component, for t = 170 
and t = 250, respectively. Solid (blue) and dashed (red) lines rep- 
resent the numerical and analytical results, the agreement between 
which is excellent. 



subsonic dark-bright soliton shortly after its generation, and 
the soliton propagates undistorted thereafter over the spatial 
region occupied by the ^-component. 

To get a deeper insight into the solitary-wave formation, we 
apply an asymptotic multiscale expansion method to Eqs. © 
and ([5]). In particular, considering only the flat region of the 
^-component (where the soliton formation is observed), we 
neglect the trapping potentials Vi } 2 in Eqs. ©-Q (recall that 
V2 is switched off in the region where V\ = and the u- 
component is practically flat), and assume that chemical po- 
tential H2 is small. We thus introduce a formal small parame- 



ter e which measures the smallness of the chemical potential 
of the v-component, and substitute /12 — » ep2 in Eqs. dU)-©. 
Our main goal is to reduce the original system of Eqs. ©-© 
to a simpler one, namely the completely integrable Mel'nikov 
system 112611 . which will provide an analytical description of 
the vector soliton observed in Fig. [4] 

The analysis starts by introducing the ansatz, 



/Mi 



qexjp(i6 



pi 



(7) 
(8) 



where real functions p(x, t) and <p(x, t) denote the amplitude 
and phase of the ^-component, while q(x, t) and 9(x, t) rep- 
resent a (complex) amplitude and phase of the v-component, 
with the real parameter Vt setting the frequency of the v- 
component. As we will see below, ft is related to the am- 
plitude of the dark soliton in the ^-component. 

Substituting Eqs. (H])-® into Eqs. (lU)-©, we obtain the 
following system of equations for fields p, (p and q: 



p t +p + P\q\ 2 -p!+p 2 x 



P~ 1/2 (P 1/2 ) XX 



0, (9) 



1 



Pt + (P¥>x)x = 0, 



i (q t + ^/2plq x ^ + q xx ~ [\q\ 2 ~ P(Pi ~ p)}Q 
+ e(fjL 2 + ^)<2 = 0. 
Next, we define the slow variables: 

T = s 3/2 t, X = s 1/2 (x-c s t), 



(10) 



(11) 



(12) 



where c s is an unknown velocity, to be determined in a self- 
consistent manner. Additionally, we express unknown fields 
p, (f and q as asymptotic series in e: 



p = p 1 +sp 1 (X,T)+s 2 p 2 (X,T)- 

- + 

q = e qi (X,T)+s 2 q2(X,T) 



ip = e 1 ' Vi (X, T) + e 3 / 2 <p 2 (X, T) + . . . , (14) 



(13) 
(14) 
(15) 
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FIG. 6. (Color online) The evolution of the largest density dip in 
the ^-component [the upper (blue) line] and peak density of the v- 
component [the lower (red) line], for the trap's velocity v = 0.84. 



Substituting expansions (TT3b-(IT4b into Eqs. d9b-dTQb. at 
the leading-order of approximation, i.e., at orders 0(e) and 
0(s 3 / 2 ), compatibility conditions of the resulting equations 
lead to the following results. First, we determine the unknown 
velocity c s = V2c s , which turns out to be the speed of sound 
of the ^-component: 

c s = ±^Ti, (16) 

and derive an equation connecting unknown phase cpi with 
amplitude pi : 



Px c s pi. 



(17) 



Proceeding to the next orders, i.e., 0(e 2 ) and 0(e 5 ^ 2 ), the 
compatibility condition for the respective equations ensuing 
from Eqs. (l9l)-([T0l) yields the following nonlinear equation: 

6 c 

2/>it H pipix - ttPixxx + c s /3(|<?i| 2 )x = 0. (18) 

c s 2 

At the order G(e 2 ), Eq. (ITTI) is reduced to the following equa- 
tion: 



qixx ~ fipiqi + O2 + tyqi = 0. 



(19) 



The system of Eqs. dT8b-dT9l) represents the Korteweg-de 
Vries (KdV) equation with a self-consistent source, which 
is determined by a stationary Schrodinger equation, is the 
Mel'nikov system [26], which is completely integrable and 
possesses exact soliton solutions of the form: 



Pi 



^sech 2 /?, 



qi — frexp (—iuJoT) sech??, 



(20) 
(21) 



where 77 = VS(X + VT), with V being the soliton velocity, 
and ujo is an arbitrary constant. Further, parameters S = P2 + 
ft and b (which set the amplitudes of pi and qi, respectively), 
are related to the soliton velocity via the equation: 



AS 



(V - Sc s ) . 



(22) 



Obviously, the velocity of the Mel'nikov's soliton is bounded 
from below, viz., V > V cv = Sc s , which follows from the 
condition that b 2 is real. 



FIG. 7. (Color online) The same as in the top panel of Fig. HI but 
for the trap's velocity v = 1. In this case, although the v-component 
is quickly dispersed (and dark-bright solitons cannot be formed), it 
creates a robust dark soliton in the ^-component. 



Using the above results, we can now write down an approx- 
imate [valid up to order O(e)] soliton solution of Eqs. ©- 
0. This solution is a vectorial soliton, composed of a dark 
soliton in the ^-component, coupled to a bright soliton in the 
v-component. In terms of the original variables x and t, it is 
expressed as follows: 

u(x,t) « ^/ii — ^sech 2 ^ exp (^i^/eSc' 1 tanhr/^ , (23) 
v(x, t) ~ ebsechr] exp(z^), (24) 

where 77 = \fs8(x — v so \t), with ^ so i = c s — eV being the ve- 
locity of the soliton. Note that this solution is characterized by 
two independent parameters, which set the dark- and bright- 
soliton's amplitudes, eS and eb. These parameters determine 
the soliton 's velocity too, via Eq. (l22b . 

This analytical result is compared to results of the numer- 
ical simulations in the middle and bottom panels of Fig. |4j 
upon numerically determining the soliton 's amplitudes, and 
finding its velocity by means of Eq. (l22l) . we can analytically 
determine the soliton density profiles [dashed (red) lines] and 
compare them to the ones found by the direct numerical inte- 
gration of Eqs. ©-(15]) [solid (blue) lines]. In this way, we find 
that the dark and bright soliton's amplitudes are s5 = 0.27 
and (eb) 2 = 0.08, respectively, and the numerically found 
soliton velocity is 0.88. The latter value is in excellent agree- 
ment with the analytical prediction, ^ so i = 0.89. This is also 
directly observed in Fig. |U where the analytical predictions 
[dashed (red) lines] are compared to results of the simula- 
tions [solid (blue) lines], see the snapshots corresponding to 
t = 170 (middle panel) and t = 250 (bottom panel). A simi- 
lar excellent agreement is observed in Fig. [5] but for a longer 
propagation time (for trap Vi of a larger size, L Xl = 123 pm, 
or L Xl = 490 in the dimensionless notation). We note in pass- 
ing that we have also checked that the Mel'nikov dark-bright 
soliton remains robust after its reflection from the trap bound- 
aries and its interaction with the sound waves. 

Coming back to Fig.O it is observed that, close to the right 
edge of the intermediate velocity regime, viz., at 0.8 < v < 
0.95, the solitons undergo a transient behavior: their peak den- 
sities perform breathing oscillations, as shown in Fig. for 
v = 0.84. Clearly, this type of the behavior cannot be cap- 
tured by the Mel'nikov model (which does not give rise to 
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breathing modes). After this transient region, in the regime 
of a large trap velocity, v > 0.95, the v-component does not 
stay localized and is quickly dispersed; this can also be ob- 
served in Fig. [2] Thus, in this case, the dark-bright soliton 
of the Mel'nikov type is not formed. However, shortly after 
penetrating the ^-component, the v-component creates a dark 
soliton there, which propagates robustly all over the spatial 
region occupied by the ^-component, cf. Fig. [7] for v = 1. 

To summarize, in this case of the v-component trap was 
switched off, we have identified three basic regimes. For 
small initial trap's speeds, the u and v components are too 
slow for combining into one of the higher- speed Mel'nikov- 
type dark-bright solitons. For the trap's speed which is too 
large, the intrusion of the v component into u produces a sin- 
gle dark soliton, but subsequently the v component disperses 
fast, without creation of other coherent structures. In the inter- 
mediate regime, the two components can "lock" into a dark- 
bright solitary wave, which can be mapped into the solitons of 
the Mel'nikov model. 



B. Evolution of the small component in the presence of the 
trap. 

We now consider the case when trap V2 is not switched- 
off and, thus, the v-component evolves in the presence of its 
trap. The pertinent contour plots depicting the largest dip in 
the ^-component and the peak density of the v-component 
are shown in the top and bottom panels of Fig. [U respec- 
tively. Similarly to the case without the trap, three different 
regimes can be identified. In particular, for low trap velocities 
v < 0.58, the dynamics of the two components is similar to 
that observed in Fig. [3] Again, in this case, long-lived local- 
ized structures are not formed in either component, and strong 
emission of radiation is observed. 

For the trap velocity in the intermediate region of 0.58 < 
v < 0.87 [see the region in between the two dashed (green) 
lines in Fig. [5), we find that dark-bright solitons of the the 
Mel'nikov type are again formed (as in the previous case 
where the trap V2 was switched off), but they eventually de- 
cay. A typical example of this behavior is shown in Fig. [51 
for v = 0.8. As is observed in the inset of the figure, after 
the v-component has intruded into the v-field, it splits into 
fragments, following its in-trap dynamics. The largest one 
among these fragments couples to a dark soliton in the u- 
component, and the resulting structure can be categorized as 
a Mel'nikov soliton: for the parameters used in Fig. [9] the nu- 
merically found average dark and bright soliton 's amplitudes 
are eS « 0.26 and (eb) 2 « 0.08, which correspond to the soli- 
ton's velocity [found via Eq. (l22h l v so \ « 1.1. This is in good 
agreement with the numerically found mean soliton 's veloc- 
ity, which is « 1, a fact that indicates that this structure is in- 
deed proximal to a Mel'nikov soliton. However, in the case of 
Fig-El me v-component travels with speed v = 0.8, while the 
soliton propagates with a larger velocity. Therefore, the con- 
fined bright- soliton component inevitably collides with - and 
is reflected by - the trap boundaries. This results in complex 
motion of the dark-bright soliton, which cannot be sustained 




FIG. 8. (Color online) The same as in Fig. [2] but in the case when 
the v-component is subject to the action of the trap. In this case too, 
three different regimes can be identified: a) < v < 0.58 with no 
localized excitations, b) 0.58 < v < 0.87, where unstable dark- 
bright solitons are observed, and c) v > 0.87, in which case we 
observe the nucleation of multiple dark solitons in the ^-component. 




FIG. 9. (Color online) The same as in the top panel of Fig. |4] but 
for the trap's velocity v = 0.8. In this case, a dark-bright soliton of 
the Mel'nikov type is formed, but it decays due to the fact that the 
soliton' s velocity is larger than the speed of the trap (see text). 



for long times and eventually decays (at t ~ 200). In fact, it 
is the trap itself in this case which is detrimental to the exis- 
tence of the Mel'nikov soliton. While the latter has its intrin- 
sic velocity upon the formation, the trap bears its own distinct 
velocity, and, barring a non-generic scenario when these two 
velocities coincide, we cannot expect the solitary wave to per- 
sist, due to the unavoidable collision with either the "front" 
or the "back" of the trap. This is the key distinguishing fea- 
ture between the intermediate speed cases in the presence and 
absence of the trap. 

When the trap's velocity exceeds the value v = 0.87, we 
observe a constant robust large dip, in the ^-component, while 
the v-component features a peak, oscillating around a rela- 
tively large value. As we show below, this is a critical ve- 
locity above which multiple dark solitons are generated in the 
^-component, indicating breakdown of the superfluidity. As 
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FIG. 10. (Color online) The top panel is the same as in Fig. [4] but 
for the trap's velocity v = 1.2. In this case, the superfluidity of the 
^-component breaks down, and multiple dark solitons are formed. 



shown in the top panel of Fig. [8l a persistent, almost zero- 
density region arises in the ^-component. In this case, the 
nucleation of multiple dark solitons in the ^-component oc- 
curs. A typical example is shown in the top panel of Fig. [TO] 
(for v — 1.2): shortly after the v-component has intruded into 
the flat part of the ^-component, dark solitons are repeatedly 
emitted. This behavior can be understood as the breakdown 
of the superfluidity in the ^-component, caused by the super- 
sonic motion of the v-component: in fact, as shown in the 
bottom panel of Fig. [TO] at early stages of the evolution, the 
v-component actually becomes a sharply localized "obstacle", 
which moves with a velocity exceeding the critical one (which 
is, roughly, the speed of sound c s of the ^-component). Ac- 
cordingly, the dynamics follows the scenario revealed in pre- 
vious theoretical J27l |28|] and experimental [Hill works con- 
cerning the breakdown of the superfluidity in quasi- ID BECs. 
In all these works, continuous emission of dark solitons was 
reported. The mechanism of the breakdown was elucidated 
in the pioneering work of Hakim ll27ll . In particular, up to 
the critical point, the moving defect (in our case, the effec- 
tive one, represented by the v-component) can sustain a pair 
of localized states co-traveling with it (a saddle and a cen- 
ter). Past the critical point, a saddle-center bifurcation arises, 
and such states can no longer be sustained. To alleviate the 
local super-criticality (in terms of the speed), the defect nu- 
cleates a (sub-sonic) dark solitary wave. Yet, once the wave 
has traveled sufficiently far downstream, it no longer renders 
the vicinity of the defect sub-critical, and a new dark soliton is 
emitted. The repetition of this process leads to the emergence 
of a dark- soliton train. 

It is relevant to quantitatively justify the above arguments, 
upon estimating the critical velocity needed for the breakdown 
of the superfluidity. This can be done as follows. First, we ob- 
serve, in the example depicted in Fig. [TT] for v = 0.8, that 
the v-component has initially (at t = 0) the form of an ex- 
tended wave packet (shaped by trap acting on it); however, 
after intruding into the v-field, it splits into fragments. The 
largest among these fragments has the peak density approx- 
imately four times as large as the initial one, and becomes 
sharply localized, with the half- width at half-maximum of the 
order of healing length £ (which is equal to 1 in our dimen- 
sionless notation). Based on this observation, it is possible 
to consider the v-component as a delta-like obstacle of some 
effective strength v* . To calculate it, we fit the v-component 
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FIG. 11. (Color online) The plot shows the density of the v- 
component for the trap's velocity v — 0.8, in two instants: at t — 
(the blue line) and at t — 32 (the green line), when it attains the max- 
imum peak density. The inset shows the density profile of the larger 
part of the v-component at t = 32; its half-width at half-maximum 
(HWHM) equals to 1.3, being on the order of the healing length, £. 



by a normalized Gaussian function with an equal amplitude at 
t ~ 32 (when the v-component attains its maximum peak den- 
sity). Then, we use the approximate expression for the critical 
velocity of Ref. [27], which reads: 



4(1 - vl/2) 



[yr+4^-(i+f c 2 f )] 1/2 

2vl - 1 + 



(25) 



In Fig. [321 we show the dependence of the critical veloc- 
ity v cr , calculated as explained above, on the trap's velocity 
v. By applying a straight-line fit to the values of v cr , we find 
that v cr « 0.048^ + 0.832. Thus, there is a very weak de- 
pendence of the effective strength of amplitude of the steep 
moving obstacle, represented by the v-component, on the ini- 
tial trap's speed. By comparing this fit with the diagonal in 
the (positive quarter-) plane of v cr vs. v, we can conclude 
what initial speeds of the trap turn out to be supercritical. In 
particular, it is clear from Fig. [12] that, for v > 0.874, the 
trap exceeds the critical velocity v CY , hence the superfluidity 
is expected to break down. Note that the numerically obtained 
velocity, beyond which the superfluidity breaks down indeed 
and dark solitons are emitted was found to be v = 0.87 [see 
the rightmost dashed (green) lines in Fig. [8), which is in ex- 
cellent agreement with the theoretical approximation based on 
Eq. ([25]). 

In summary, three regimes are also identified in the case 
when the v-component is subject to the action of the poten- 
tial, although their dynamics is different from the case when 
the trap was absent, except for the low-speed case. In the 
intermediate- speed regime, solitary waves of the Mel'nikov 
type form but are unstable due to the velocity mismatch with 
the moving trap. In the large-speed regime, the superfluidity 
was shown to break down at the critical point very close to 
the related analytical prediction, which leads, as in the case 
of Ref. ll27ll . to the emission of a train of dark solitons, due 
to the effective potential exerted by the second (v) component 
on the first (u) one. 
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FIG. 12. (Color online) The critical velocity, v C v, versus the trap's 
velocity, v. First we calculate the value of v cr [depicted by (blue) 
circles] for different trap's velocities, and then apply the straight-line 
fit [the dashed (green) line]. Breakup of the superfluidity is expected 
when the trap's velocity becomes larger than the critical velocity, i.e. 
at v w 0.874, as indicated by the vertical dashed (black) line. 



IV. CONCLUSIONS. 

In this work, we have studied the dynamics of two quasi- 
1D counter-propagating immiscible superfluids. Our setup is 
based on the BEC mixture composed by two different spin 
states of the same atomic species, with the number of atoms 
in one of the two components being much smaller than in the 
other. In our simulations of the coupled GPEs, we have as- 
sumed that the "large" component is at rest while the "small" 
component starts intruding into the larger one, moving at a 
constant velocity v. 

We have considered two different situations, namely, with 
the small component propagating either in the absence or in 
the presence of the trap. For sufficiently small trap's velocities 
v, we have found that no localized structures emerge. 

On the other hand, for intermediate velocities in the in- 
terval of 0.61 c s < v < 0.95c s , if the trap acting on the 
small component is switched off, we have found that a robust 
dark-bright soliton is formed. This structure is composed of a 
bright (dark) soliton in the small (large) component, and can 
be very well approximated analytically by means of the mul- 
tiscale asymptotic expansion method, which reduces the orig- 
inal GPE system to the completely integrable Mel'nikov sys- 
tem. This system has exact soliton solutions that can be used 
to construct the approximate dark-bright soliton solutions of 
the GPE system. We have also found that, if the small compo- 
nent is subject to the action of the trapping potential, then the 
Mel'nikov-type dark-bright soliton forms in a similar velocity 
regime, 0.58c s < v < 0.87c s , but it cannot be sustained due 



to the mismatch between its own velocity and the speed of the 
moving trap. 

If the small component propagates in the absence of the 
trap, then, for large initial trap velocities, v > 0.95c s , it dis- 
perses and no dark-bright soliton is formed; nevertheless, a 
dark soliton is always created in the large component. On the 
other hand, if the small component propagates in the presence 
of its trap, which moves with speed v > 0.87c s , we have 
found that it gets deformed into a sharply localized object of 
the width on the order of the healing length, which propagates 
with a velocity larger than the critical velocity of the large 
component. Approximating the small component by an effec- 
tive delta-like obstacle acting on the large component, we have 
found that, in the case of such large trap velocities, the dynam- 
ics of the system follows the scenario known for the quasi- ID 
superfluid flow past an obstacle. In this way, the superfluidity 
of the large component breaks down, and nucleation of a train 
of dark solitons occurs. 

Our results suggest that such counterflow settings may be 
quite relevant to the observation of fundamental phenomena, 
such as the formation of solitons and/or the identifying the 
critical velocity in superfluid flows. It would be interesting to 
extend the considerations to higher-dimensional settings and, 
also, to larger number of components, as, e.g., in the case of 
F = 1 or F = 2 spinor BECs 113 311 . In that case, vortex- 
bright solitary waves [l34ll and their dipole-more generaliza- 
tions II 3 511 are expected to arise, as well as their spinor coun- 
terparts. The generalizations with a larger number of com- 
ponents and in higher dimensions is likely to produce novel 
waveforms. Another relevant extension concerns the role of 
gi2 (the cross-component interaction strength) and the ap- 
proach to the miscibility-immiscibility threshold. In that con- 
text, exploring the conditions under which additional wave- 
forms (such as dark-dark solitons of [[TBI [l7|] and their multi- 
dimensional counterparts) can arise would also be a theme of 
interest. These topics are currently under investigation and 
will be reported elsewhere. 
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